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Theory of magnetic deflagration (avalanches) in crystals of molecular magnets has been developed. 
The phenomenon resembles the burning of a chemical substance, with the Zeeman energy playing 
the role of the chemical energy. Non-destructive reversible character of magnetic deflagration, as 
well as the possibility to continuously tune the flammability of the crystal by changing the magnetic 
field, makes molecular magnets an attractive toy system for a detailed study of the burning process. 
Besides simplicity, new features, as compared to the chemical burning, include possibility of quantum 
decay of metastable spin states and strong temperature dependence of the heat capacity and thermal 
conductivity. We obtain analytical and numerical solutions for criteria of the ignition of magnetic 
deflagration, and compute the ignition rate and the speed of the developed deflagration front. 

PACS numbers: 75.50.Xx, 76.60.Es, 82.33.Vx 
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I. INTRODUCTION 

Recently, it has been observed that molecular mag- 
nets exhibit explosive relaxation towards thermal equi- 
librium that resembles propagation of a flame through a 
flammable chemical substanceJ^ Theory of this effect is 
the subject of this paper. 

Crystals of molecular magnets first attracted attention 
of physicists after it was demonstratedii that individ- 
ual molecules inside such crystals behave as superpara- 
magnetic particlesji^ Due to large molecular spin (e.g., 
S = 10 for Mn-12 and Fe-8 molecular magnets) and 
high magnetic anisotropy, spin-up and spin-down states 
of many molecular magnets are separated by a large en- 
ergy barrier. Consequently, unlike conventional param- 
agnets, molecular magnets are characterized by a macro- 
scopic time of thermal relaxation between spin-up and 
spin-down states. Similarly large times are needed for 
quantum transitions between these states to occur, which 
allows one to speak about quantum tunneling of the mag- 
netic moment.^ Due to this effect, molecular magnets ex- 
hibit spectacular staircase magnetization curve?^ It has 
been known for some time that the low-temperature mag- 
netic relaxation in crystals of molecular magnets can oc- 
cur via two mechanisms. The first, slow mechanism in- 
volves random thermal and quantum spin transitions at 
unrelated spatial points. Such transitions influence each 
other only through weak long-range dipolar fields asso- 
ciated with the magnetic moments of the moleculeSfi^ 
The resulting relaxation lasts macroscopic times. This 
allows one to study transition rates by simply measuring 
the time dependence of the macroscopic magnetization 
of the crystal. 

The second mechanism of relaxation - magnetic 
avalanches - corresponds to the abrupt reversal of 
the magnetization when a sufficiently large crystal is 
placed in a large magnetic field opposite to its magnetic 
moment i^iSiiii The avalanche was long believed to be a 
thermal runaway in which the Zeeman energy released 
by the relaxing molecules gets transformed into heat that 
generates transitions in the neighboring molecules and 



accelerates the total energy release. Such a relaxation 
that typically occurs in a millisecond time was long con- 
sidered a nuisance as it often interfered with experimen- 
tal studies of spin tunneling. More recently, it was re- 
alized through time-resolved local measurements^^ that 
magnetic avalanches resemble propagation of a flame - 
deflagration - in which the role of the chemical energy 
stored in a flammable substance is played by the Zeeman 
energy. Due to quantum tunneling between spin states 
that occurs at discrete values of the magnetic field, mag- 
netic defiagration also exhibits quantum features.-^ 

Experiments performed to date have established with 
certainty that magnetic avalanches in crystals of molec- 
ular magnets correspond to the propagation of a nar- 
row front of the magnetization reversal. The analogy 
with burning of a flammable chemical substance has been 
confirmed by the study of the dependence of the flame 
speed on the energy barrier. In crystals of molecular 
magnets, the latter can be continuously tuned by the 
magnetic field. For the study of deflagration this tun- 
ability of the barrier, as well as the reversible nature of 
the magnetic burning, provides a great advantage over 
irreversible burning of a chemical substance with a fixed 
energy barrier. Thus a detailed study of magnetic defia- 
gration can answer important questions of the theory of 
combustion and detonation. ^'^ There are also novel fea- 
tures that are absent in conventional combustion. They 
include a very strong temperature dependence of the spe- 
cific heat and thermal conductivity of molecular magnets 
at low temperature-''^'-^ and the possibility of magneti- 
zation reversal via quantum tunneling. 

In this paper we intend to answer the following ques- 
tions: 

• The critical combination of parameters (magnetic 
field, initial temperature, and the size of the sam- 
ple) that sets off the deflagration process. 

• The mode of instability. 

• The time that elapses between bringing the system 
above the deflagration threshold and the ignition of 
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FIG. 1; Energy of a molecular magnet as function of Sz. 
Quantum energy levels are shown by black circles for 5* = 4. 



the deflagration process (the ignition time). 

• The temperature of the flame and the velocity of 
the deflagration front. 

We will show that the ignition of magnetic deflagration 
in molecular magnets is very different from the ignition 
of magnetization reversal in ferromagnets. The latter is 
dominated by the exchange interaction and begins with 
the nucleation of a small critical nucleus of opposite mag- 
netization that spreads and occupies the entire sample. 
On the contrary, the magnetic deflagration in a param- 
agnetic crystal of magnetic molecules begins as a large- 
scale instability of a smooth temperature profile inside 
the sample against formation of a rapidly moving defla- 
gration front. 

The structure of the paper is as follows. Properties 
of molecular magnets will be discussed in Sec. |TT1 The 
mechanism of thermal runaway in a crystal of magnetic 
molecules will be analyzed in Sec. IIIII Stability of the 
quasi-stationary temperature profile in a crystal of molec- 
ular magnets will be analyzed in Sec. IIVI The ignition 
rate will be studied in Sec. |Vl Structure and the velocity 
of a developed deflagration front will be investigated in 
Sec. I VII Numerical illustrations of the deflagration pro- 
cess will be given in Sec. IVIII Relevance of our results to 
experiment and possible future directions of theory and 
experiment will be discussed in Sec. IVIIII 



II. MOLECULAR MAGNETS 

A. Magnetic bistability and spin tunneling 

A single molecule of a molecular magnet can be de- 
scribed by the Hamiltonian 



n = -DSl-giiBH,S,+n\ 



(1) 



where S is spin, D is the constant of the uniaxial z 
anisotropy that creates magnetic bistability, is the 



bias magnetic field, and Ti.' is a small part of the Hamil- 
tonian that does not commute with Sz and is responsi- 
ble for spin tunneling. If S is large (as, e.g., in Mni2 
and Fcg), magnetic bistability can be reasonably well de- 
scribed within the classical model with the energy that 
depends on the classical vector s = S/S" and has the form 



Here 



2hsz) Uo. 



_ gusHz 



2DS 



(2) 



(3) 



are the zero-bias energy barrier and the reduced bias 
field. The dependence E{sz) is shown in Fig. [1] The 
spin-projection value corresponding to the barrier be- 
tween the two wells follows from dE/dsz = and is given 



by s 



-h. The minima of E and its value at the top 



of the barrier are 



E± = -{l±2h)Uo, Eb^h^Uo. 



(4) 



Thus the values of the energy barriers for the molecules 
on the left and on the right are given by U± = Eh — E± = 
(1 ± hy Uq. In the case of /i > 0, that we will consider 
throughout the paper, Sz = — 1 is a metastable minimum, 
whereas 1 is the absolute minimum of the energy. 

Below we will use U- = U, 



U={1- hf Uo. 



(5) 



The energy difference between the two minima is given 
by 



AE = E^ - E+ = AhUo. 



(6) 



The noncommuting term TC in Eq. ([T]) gives rise to 
resonance spin tunneling between the states at the two 
sides of the barrier if the bias field satisfies the condition 



g^lBHz = kD, 



0,±1,±2, 



(7) 



This modifies the process of thermal activation of spins 
at low temperatures. Off resonance, the spins have to be 
thermally activated all the way up from the bottom of 
the metastable well to the top of the barrier. On reso- 
nance, however, it is sufficient to be thermally activated 
up to the energy level below the barrier where spin tun- 
neling is sufficiently strong to take the molecule to the 
other side of the barrier. This leads to the resonance in- 
crease of the relaxation rate, see Fig. 7 of Ref. ^ On the 
phenomenological level, this effect can be encapsulated 
into the effective barrier with dips at the resonance bias 
fields given by Eq. ([7]). Since the exact form of the ef- 
fective barrier depends on the form of Ti' that we do not 
analyze in this paper, we will use for numerical work the 
fitting function taken from experiment a^'^^ on Mni2 and 
replace Eq. ^ by U{h) = u{h)Uo with 



u{h) = (1 - /i)- 0.0806 



1 



^ D 



(8) 
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B. Magnetic relaxation and heat transfer 



At low temperatures 



U 



VF> 1 



(9) 



thermally activated transition of magnetic molecules over 
the barrier can be described by the rate equations for the 
numbers of molecules in the left and right wells, n±, that 



satisfy + n_ 



1. The equation for the number of 



particles in the metastable well n_ has the form 



,(eq) 



(10) 



where F — + F ^. In accordance with the detailed- 
balance condition 



,(cq) 
^ — h n_ 



cxp 



AE 



Acq) 



1 



exp 



(11) 



Using F-i = Foe for W ^ 1, one obtains 



r = roe 



1 + cxp 



AE 



(12) 



where the second term square brackets describes transi- 
tions from the stable well to the metastable well. In the 
strong-bias case, AE ^ ksT, this term can be omitted. 
This yields simply F = Foe"^, and sets nL^'^^ = (full 
burning) . 

When a magnetic molecule makes a transition from 
the metastable state Sz = —1 to the absolute energy 
minimum = 1, the energy AE is released. Thermal- 
ization of this energy leads to the temperature change 
AT — AE/Cph, where Cph is the phonon heat capac- 
ity per magnetic molecule. Other contributions to the 
specific heat at low temperatures are considered small. 
The magnetic relaxation creates a source in the heat con- 
duction equation. Another term in this equation is the 
divergence of the heat flow 



(13) 



where k is thermal conductivity. The full system of equa- 
tions for the temperature T and the population of the 
metastable minimum n_ has the form 



dT 
'dt 



1 

Cph 
-F 



V • fcVT 



AE dn- 
Cph dt 



Am) 



(14) 



An important feature of magnetic deflagration is 
strong temperature dependence of the heat capacity and 
thermal conductivity at low temperatures. As the tem- 
peratures before and behind the deflagration front can 
differ by an order of magnitude, this effect cannot be 
neglected. The phonon heat capacity Cph has the form 



Cph Aki 



T 



(15) 



where a = 3 in three dimensions, A is a numerical factor 
and Od is the Debye temperature. At low temperatures 
only acoustic phonons are excited, whereas high-energy 
optical phonons are frozen out. Thus one can use the A 
value for the simple model of a crystal,'* A — 127r'*/5 ~ 
234 that is in a good agreement with measurements- on 
Mni2. The thermal diffusivity 



K — k/Cph 



(16) 



depends on the average mean free path of thermal 
phonons. At low temperatures the main scattering mech- 
anism is scattering on impurities, so that (see Ref. i7| and 
references therein) 

KcxT^^, /3 = 13/3. (17) 

Accordingly the thermal conductivity behaves as 



kcxT- 



7 = /5- 



4/3 



(18) 



for a = 3. 

The heat-conduction equation can be brought into a 
more elegant form by choosing the phonon energy £ as 
the dynamical variable. Using Cph = d£ /dT one obtains 



'dt 



V • kSIE - AE 



9n_ 



(19) 



This form of the equations is convenient for the study of 
the stationary deflagration front as it allows one to im- 
mediately obtain the first intergal of the heat-conduction 
equation. 

Alternatively one can use 



K = I k{T')dT' 

J To 



(20) 



as the temperature variable, with Tq being a reference 
temperature. With this choice, the first of Eqs. takes 
the form 



IdK o , an_ 

= V K - AE , 

K dt dt 



(21) 



while in the second equation one should use T — T{K) in 
the expression for F. This form of equations is convenient 
for the study of the deflagration threshold in the case 
when the temperature along the boundary of the crystal 
varies. 



III. THERMAL RUNAWAY 

Deflagration begins with a thermal runaway in a part 
of the sample that has a lower barrier [/ or a higher tem- 
perature T than the surrounding area. Thermal runaway 
needs some time to develop. We call it the ignition time 
Tig. The shortest ignition time is achieved if the heat re- 
leased by the relaxation remains in the sample and does 
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not escape through its boundaries. We denote this igni- 
tion time tI^\ as it is related to the ignituin in the in- 
finite sample, see below. If the rate of heat transfer out 
of the sample is sufficiently high, the ignition does not 
occur. This explains why small crystals do not exhibit 
magnetic avalanches. In this section we derive the expres- 
sion for the ignition rate rl^-* = l/r-^^ that plays fun- 
damental role in subsequent considerations and justifies 
the validity of the explosive approximation n_ =^ n__i 
(n_,i is the initial value of n_) that will be used below. 

For an infinite and/or thermally insulated sample one 
can drop the diffusion term in the first of Eqs. p^ . This 
yields 



dT _ AE 



r{T)n-. 



(22) 



The initial conditions are T = Tq and n_ = ri-^i. Using 
Eq. in the strong-bias case, it is convenient to intro- 
duce the reduced temperature deviation and the reduced 
population of the metastable well 



Wo- 



T-Tn 



To 



(23) 



with Wq defined by Eq. ([9]) with T — Tq. Linearization of 
the argument of T{T) on 6 leads to the system of equa- 
tions 



89 



dn 



— = e n, 1^0 — 

OT OT 



(24) 



is given by 



where r = t is the reduced time. The ignition rate 



where 



1^0 = Wo 



Cph,o2o 



u a 



ph,0 



(25) 



(26) 



and Cph.o is the phonon heat capacity at T = Tq. Note 
that in cases of practical interest vq is a large parameter 
since it contains Wq, whereas AE/U is typically of order 
one. The factor fcB/Cph,o is also large at small temper- 
ature To, see Eq. (ITS]) . Since vq 1, the evolution of 
n_ is much slower than that of 6, so that one can ap- 
proximately replace n_ with to study ignition. This 
is what we call the explosive approximation or the ex- 
plosive limit. (Indeed, for explosives the parameter vq is 
very large, so that the problem of their stability can be 
considered without taking into account that a small part 
of the explosive has already burned and the heat release 
has been reduced because of this.) 

Setting n_ in the first of Eqs. one ob- 

tains an isolated equation that has the solution 9{t) = 
— ln(l — r) reaching infinity exactly at r = 1. In real 

units, it corresponds to t — t|^^ Of course, soon after 
the ignition the deviation of T from Tq becomes large, 
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FIG. 2: The dependence of the magnetic field that ignites 
deflagration on the temperature of the crystal of Mni2. Di- 
mensionless parameters are given by Eqs. (|3]) and pO|l. 



thus the linearized equation drO — , as well as the re- 
placement n_ =5> Ti^_i, becomes invalid. Nevertheless, 
the dominant contribution into the ignition time comes 
from the time range 0^1 where still T — Tg <C To and 
the equation drO = is valid. 

The low-tcmpcraturc explosive approximation n_ 
n^_i introduced above drastically simplifies the problem 
of the ignition of deflagration and allows one to under- 
stand it in simple terms. With ri_ = n^j, the tempera- 
ture of the sample is the only relevant variable, and its 
dynamics is determined by the competition of the two 
terms. One of them is the heat release due to the re- 
laxation that is strongly nonlinear in temperature. The 
other one is the heat loss due to the heat conduction that 
is linear on temperature but contains spatial derivatives. 
There are two scenarious if one starts with the sample 
having an uniform temperature T = Tq that coincides 
with the constant temperature of the sample boundaries. 
In the first scenario, the temperature increases because 
of the heat release, typically with a maximum at the cen- 
ter of the sample, until a sufficient temperature gradient 
develops that provides the balance between the heat re- 
lease and heat loss through the boundaries. The resulting 
state is the stationary state of the system. In the second 
scenario, the heat loss through the boundaries is insuffi- 
cient to balance the increase of the heat-release due to the 
rise of temperature. This happens, in particular, if the 
sample is sufficiently large. In this case there is no sta- 
tionary state, and the temperature growth, slow at the 
beginning, leads to a thermal runaway. Changing one 
of the parameters (sample size, energy barrier, temper- 
ature at the boundaries, initial magnetization) one can 
reach the situation in which the stationary state disap- 
pears and the runaway begins. We call it the ignition 
threshold. Below we present analytical and numerical 
results for the ignition threshold in different cases in the 
explosive limit. We will see that the ignition is mainly 



5 



controlled by the parameter 



(27) 



where i? is a typical shortest distance from the center of 
the sample to its boundary and fco is thermal conductivity 
A; at T = Tq. For a Id sample (i.e., a cylinder thermally 
insulated along its side) one has R = L/2, where L is the 
sample thickness, whereas for cylindrical and spherical 
samples R is the radius. r|°°^ is given by Eq. (p5)) while 
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Finding 8^, for t/ = 1 , 2, 3 















2ko/R^ 



(28) 



is the rate of thermal equilibration within the sample at 
T — Tq.Ixv Eq. (P7|) Iq is the characteristic thermal length 



(oo) 



(29) 



at T = Tr 



0- 



We will see that in the simplest case of the uniform en- 
ergy barrier and constant temperature Tq maintained at 
the boundaries (uniform conditions) the ignition thresh- 
old corresponds to (S = (5c ~ 1. The exact value 5c de- 
pends on the geometry of the sample. In particular, in 
one dimension 5c — 0.439. This allows one to obtain a re- 
lation between the temperature Tg and the barrier U at 
the ignition threshold, that also depends on the sample 
size R and other parameters. In terms of dimcnsionless 
parameters 



A: 



R^Tpn-,, 
2k{fo)/ki 



(30) 



[see Eqs. ^ and p^ ] one can write Eq. (P7)) in the form 



u{h) 



(31) 



where u{h) is given by u{h) = (1 — h)"^ for the classical 
model and by Eq. ^ with account of spin tunneling. Re- 
solving the threshold equation 5 = 5c requires the knowl- 
edge of the temperature dependence of thermal conduc- 
tivity k. It turns out, however (see below) that the exact 
form of fc(To) given by Eq. (|18p is not essential. The re- 
sults for k{T) = const and thus A = const are shown in 
Fig. [5] for three different values of A. One can see that 
increasing A (say, due to the increasing of the sample 
size R) leads to the decrease of the critical values of Tq 
and h. In the realistic case of the large Arrhenius expo- 
nent in Eq. (I3ip . the dependence on A and thus on k is 
logarithmic. 



FIG. 3: Plot of (5(Smax) for d = 1, 2, 3 that allows one to obtain 
the ignition threshold from the maximuin of these curves. 



IV. IGNITION THRESHOLD IN MOLECULAR 
MAGNETS 

A. Ignition threshold under uniform conditions 

1. Ignition threshold in one dimension 

Consider a slab of thickness L = 2R. We will see that 
at ^ 1, if the constant temperature Tq is maintained 
at the boundaries of the slab, the solution for T at the 
ignition threshold only slightly deviates from Tq. In this 
case one can linearize the problem around Tq using the 
temperature deviation 9 of Eq. ([^5]) and set k ^ kg. 
With dT/dt = in the stationary case one obtains the 
equation 



+ 25e^ = 0, 



0(±1) = O, 



(32) 



where the reduced space variable x is normalized by R. 
The first integral of this equation is 



(33) 



where 0max is the integration constant that equals to the 
maximal value of 9 achieved in the middle of the sample. 
Integrating Eq. ([33]) one obtains 



- 2 In cosh y/Se^^ 



(34) 



where the value of 0max follows from the boundary con- 
ditions 9{±.l) — 0. To find the ignition threshold, one 
can solve this equation for 5: 



In' 



1 



(35) 



The dependence 5(6'niax) is shown in Fig. [3] It has a 
maximum at 0max = ^max,c = 1.18684. The maximal 
value of 5 



5r = 0.439229, 



(36) 



6 



corresponds to the ignition threshold. Indeed, for 5 < 6c 
there are two solutions for ^maxi and the smallest of the 
two corresponds to the stationary solution of the heat- 
conduction equation. For 6 > 6c the stationary solution 
disappears. 



2. Ignition threshold in two and three dimensions 

For a cylindrical {d = 2) and spherical {d = 3) samples 
the generalization of Eq. (j32p is 



dr^ 



d-ide 

r dr 



+ 26e^ = 



(37) 



with r normalized by R and with the boundary conditions 
e'{0) = and 6{1) = 0. For d = 2 the exact solution of 
Eq. jSIl) is 



(38) 



0{r) = 21n 



1 + VT~6 + (1 - VT~6) r2 ' 



Its maximal value 



= 21n 



1 + VT~6 



(39) 



is achieved at r = 0. The ignition threshold can be found 
by the same method as in Id. Resolving this equation for 
6 one obtains 



6 = ie-'-- {^e'--/^ - 1 



(40) 



This function has a maximum at ^max — ^max ,c — 
2 In 2 ~ 1. 386 3, and the corresponding critical value of 6 
is 



6c = 1. 



(41) 



For d = 3 we are unable to find the solution of Eq. ([57]) 
in terms of known functions. Numerical solution for the 
ignition threshold in 3d consists of the following steps: (i) 
One solves Eq. (|37p with the boundary conditions 6' (0) = 
and 9{Q) — 6'max and (5 as a free parameter; (ii) One 
finds (5 as a function of ^max from the boundary condition 
6{1) = 0; (iii) One finds critical parameters from the 
maximum of (5(0max)- Our results for d = 1,2,3 are listed 
below 



(42) 



One can see that approximately 6c oc d. The curves 
6{0max.) are plotted in Fig. [3] 



B. Ignition threshold in the presence of field 
gradient 

Consider a one-dimensional problem of ignition with 
temperature at both ends maintained at Tq and the bar- 
rier U varying in space due to the gradient of the bias 



d 


6c 


6c/d 


l/6c 


^inax ,c 


1 


0.4392 


0.4392 


2.277 


1.187 


2 


1 


0.5 


1 


1.386 


3 


1.661 


0.5537 


0.6020 


1.607 



1.5 

G 



Temperature profile at criticality 
with the bias-field gradient 

P = 0, 1, 3, 10 




Heat-release profile at criticality 


g=exp[-p(l+:t)+e] 


with the bias-field gradient 






P = 0, 1, 3, 10 




o\ 






\ V\ 

io\ 









FIG. 4: Numerical results for temperature and heat-release 
profiles. 
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FIG. 5: Ignition threshold in the model with the bias-field 
gradient. 



field. Although the relative variation of U is small, the 
effect can be large for large Atthcnius factors W as its 
variation 6W can be large. Assuming that the barrier 
is the lowest at the left end of the sample and the field 
gradient is small and constant, one can write 



6W'^w{\ + x) 



(43) 
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The equation for the stationary temperature profile be- 
comes 



^ + 2^e-"'(i+-)+^ = 0, 9{±1) = 0. 



(44) 



For w ~ 1, Eq. (|32p can only be solved numerically. Here, 
instead of plotting 6 vs 6'inax, it is more convenient to plot 
5 vs 9'{—l). Numerical results for temperature profile 
9{x) and heat-release profile q{x) oc exp[— + x) + 9] 
at the ignition threshold are shown in Fig. [4l As the 
field gradient goes up, the maxima of these curves shift 
towards the end of the sample where the barrier is lower. 
The threshold condition Sc{w) is shown in Fig. \5\ The 
value of 6c increases with w since the favorable condition 
for burning is realized in a more and more narrow region 
at the left end of the sample, and there is an increasing 
heat flow out of this region in both directions. 

In the case of w ^ 1 the ignition occurs very close to 
the left end, x ~ ^1, and the heat release proportional 
to exp [—w{l + x) + 9] is very close to zero except in the 
vicinity of the left end. The temperature profile for w 3> 
1 consists of two regions: Very close to the left end the 
temperature rises sharply to the maximal temperature 
^max and then goes linearly down to zero at the right 
end. Thus for w 3> 1 one can introduce a new variable 
u = w{l + x), a new function § — 9 — u, and replace Eq. 
da by 



(45) 



with the boundary conditions i?(0) = 0, {oo) — —1. 
This equation is similar to Eq. ([5^ and its solution reads 



^{u) = l?,„ax - InCOsh^ (V^e''»- {u - M,nax)) , (46) 

where i?max and Umax are integration constants. From 
the boundary condition 'i?'(oo) = —1 one obtains 

: 1 and thus 



2V^ 



mi. 

4(5 



(47) 



Then the other boundary condition, = 0, gives 

1 



= 2 arccosh- 



Since i^max ^ 0, the ignition threshold is defined by 



(48) 



1 

4' 



w 

T 



(49) 



For 5 = Sc one has Umax = , that is, the ignition occurs 
at the left boundary. Corrections to Eq. move the 
maximal-temperature point a little away from the edge of 
the sample. Comparing this situation with the uniform 
case with the same barrier as the barrier at the left end 
(the minimal barrier), one can see that in the case of 



1/7^1 the ignition requires a much higher temperature 
To- 

The constant w in Eq. (j43| can be expressed via the 
gradient of the bias field with the help of Eq. (O: 



w=^^Wo{l~h)\Ah\^ 
dx 



(50) 



where Aft, is the change of h across the sample. In terms 
of the energy bias AE defined by Eq. ([6]) one has 



w 



1-hA [AE) 



(51) 



where A (AE) is the change of the energy bias across 
the sample. Note that at a very low temperature Tq, the 
condition w 3> 1 does not necessarily require a large field 
gradient. Because of the field gradient, the Arrhenius 
factor W increases by 1 at the characteristic distance 



^ _2R _ L _ 



L 



Wo{l^h)\Ah\ 



(52) 



from the end. It is the width of the ignition region near 
the low-barrier end. In the case of w ^ 1, that is Ih ^ L, 
the far end with the highest barrier becomes irrelevant 
for the ignition, and the ignition threshold S ^ Sc that 
follows from Eqs. (P7|) and becomes 

(53) 



Ih — io- 



For a cylinder of radius R with the bias field linearly 
changing along its symmetry axis, this condition holds if 
Ih ^ R, so that the heat flows along the cylinder axis 
z away from the face with the lowest barrier rather than 
towards the side walls of the cylinder. The problem then 
becomes one dimensional. 



C. Ignition threshold in the presence of 
temperature gradient 

While it is experimentally difficult to create a large 
gradient of the bias-field over the length of a small crys- 
tal, it is relatively easy to break the symmetry of the 
system by creating a large temperature gradient. This 
can be done by, e.g., maintaining temperature Tq at the 
left end and having Ti < Tq at the right end of a Id 
sample. If these two temperatures differ essentially, one 
has to take into account the temperature dependence of 
thermal conductivity k{T) that is strong at low T, see 
Eq. (|18p. In this case it is more convenient to use K 
defined by Eq. (|20|) instead of T. As the ignition occurs 
closer to the hot left end, it is convenient to choose Tq as 
the reference temperature and introduce 



Wn 



K 
' koTo 



(54) 



that generalizes Eq. with fco = k(To). The relax- 

ation rate can be expanded similarly to the above, 

r(r(i^)) - r(ro) exp (wo^i^) = r{To)e'. (55) 



8 



Ignition threshold for the Id problem 
with different temperatures at the ends 




0.0 



Location of the maximal temperature point 
at the ignition threshold for the Id problem 
with different temperatures at the ends 




10 15 20 25 |e I 30 



FIG. 6: Numerical results for Sc (a) and a;niax,c (b) vs the 
temperature-bias parameter Oi for the Id model with different 
temperatures at the ends, x — corresponds to the center 
of the sample. The dashed line in (a) is the asymptote 5c = 
(i6lil+21n2)Vl6 + l/4at j^ij > 1. 



The applicability of this expansion requires |<5r|/To ^ 1. 
However, in the case of Wq ^ 1 this expansion prac- 
tically works in the whole range of i5T < since the 
burning rate F becomes neghgibly small long before the 
condition \ST\/To ^ 1 is violated. In the stationary state 
9{x) satisfies the same Eq. ([5^ but with the boundary 
conditions 

0(-l) = O, 0{l) = e, = Wo^^^ <0. (56) 

In the absense of the heat release due to burning the 
solution for 9 would be a linear function, 0(x) = (1 -t- 
x)9i/2. An estimate for the width It of the region near 
the hot end where the ignition occurs can be obtained by 
setting 9{x) ~ — 1. In real units ignition occurs at the 
distance of order 



It = 4.R/\9i\ 



(57) 



where the numerical factor 4 anticipates Eq. . For the 
very cold right end, 16*11 3> 1, one has It <C R. The first 
integral of Eq. (|23p is Eq. (j33p in which the maximum 



can be achieved at some Xmax shifted from the 



central point a; = 0. The solution for 9{x) reads 

9{x) = 9max - 2 In cosh Vje^^ {x - Xmax) (58) 

which is the generalization Eq. (j34[) . Elimination of Xmax 
using the boundary conditions and solving for S yields 



6 = 



(59) 



c.f. Eq. (|35p . Computing the maximum of this function 
on ^max allows to determine Sc for any value of 9i. 

If the temperature at the right end is low, —9i is a 
large positive value, so that Eq. ([59)1 simplifies. In this 
case at the ignition threshold one has ^max <C 1, so that 
Eq. ([55)1 becomes 



1 



41n2 4^/9^, 



(60) 



The maximum of the rhs is attained at ^max — ^ma? 
4/6*1 ^ 1- Thus the ignition threshold is defined by 

'\9i\ -f 21n2 



V 4 fcoTo 



(61) 



The maximal-temperature point at the ignition threshold 
is 



= -1 



(62) 



which is close to the left end. One can see that for 16*1 1 ^ 
1 the value of 6c is large, so that much larger L is needed 
to reach the threshold for the same temperature To, as 
compared to the uniform case. The ignition threshold 
i5 = 5c is equivalent to 



It — lr\ 



(63) 



where It is given by Eq. ([57]) . This result is similar to 
Eq. ([55)1 . The numerically obtained dependences (5c(|6'i|) 
and Xmax, 0(16*1 1) are shown in Fig.[Sl 

The remaining task is to relate the temperature-bias 
parameter 9i to the temperatures at the ends, Tq and Ti. 
At a constant thermal conductivity k one obtains 



(64) 



that tends to —Wq in the limit Ti 0. Since Wq » 
1, there is a strong suppression of the ignition by the 
cold end. The effect is even stronger for the power-law 
dependence of k given by Eq. Parametrization 



7- \ -7 



(65) 



and integration in Eq. (|20|1. with the lower limit being 
To, yields 



K{T) 



koTo 



7-1 



1 



(66) 
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Now one obtains 



7-1 



(67) 




7-1 



(68) 



In the realistic case given by Eq. (fT8|) one has 7 — 1 = 1/3. 
That is, if Ti goes to zero, 9i becomes infinite negative 
and Sc diverges. This means that for a sample of any size 
one can suppress the ignition by making the temperature 
of the cold end very close to zero. This is a consequence 
of the divergence of thermal conductivity at T = 0. 



V. RATE OF IGNITION OF DEFLAGRATION 
IN MOLECULAR MAGNETS 

Above the ignition threshold, 6 > 6c, the heat loss 
via heat conduction cannot compensate the heat release 
due to burning and there is no stationary solution for 
the temperature. The temperature growth leads to a 
thermal runaway after the ignition time rjg, followed 
by the deflagration. At first we investigate the ignition 
time within the explosive approximation, n_ ^ n^ i, in 
terms of 9 defined by Eqs. (gS]) or Then we study 

deviations from the explosive approximation using a 
more general system of equations containing both 9 and 
n,-. The initial condition in all examples considered in 
this section is thermal equilibrium reached in the ab- 
sence of relaxation, for instance, a uniform temperature 
throughout the sample. This is the most transparent 
case theoretically but it may be difficult to realize 
in experiment if the ignition rate is large. At the end 
of this section we discuss other kinds of initial conditions. 

Please, get the full text of the paper here: 

www.lehman.edu/faculty/dgaranin/deflagration.pdf 



VI. STRUCTURE AND VELOCITY OF THE 
DEFLAGRATION FRONT 

A. Thermodynamics of magnetic deflagration 

Ignition of deflagration leads to a strong increase of 
the temperature and relaxation rate that results in equi- 
libration of energy between spin and phonon subsystems. 
Since deflagration is a fast process, one can use energy 
conservation, neglecting excited states and the heat loss 
through the boundaries during deflagration, 



(69) 



Here £ij = £{Tij) are the phonon energies at the initial 
and final temperatures (before the defiagration front) 
and Tf (behind the front), n_^i is the initial population 

of the metastable well, and n^'^^(T/) is the equifibrium 
value of n_ at Tf given by Eq. (fTT|) . We will call Tf the 
flame temperature. Eq. (|69p is a transcedental equation 
for the flame temperature Tf. Since Tf ^ Ti, the initial 
phonon energy £i can be neglected. Tf can be found 
analytically if £ has the form 



£ 



a + 1 [Od 



(70) 



that follows from Eq. ^ if n'il'^^Tf) is negligibly small 
— the full- burning case. In this case from £{Tf) — 
n_ , A_B one obtains 



Tf^Qi 



(a -I- 1) n^^iA.E\ 



l/(a + l) 



AkB^D 



(71) 



The full-burning condition is fcsT/ ^ AE. Using Eq. ^ 
one can rewrite it in the form of the strong-bias condition 



h > hfbn^Jj, 
where the full-burning field 



'■fb 



ks^D f a + 1 
4% I A 



l/c 



(72) 



(73) 



is a material parameter. With a = 3, A ~ 234, C/g ~ 65 
K and O/j = 38 K for Mni2 one obtains 



hfh ~ 0.038. 



(74) 



We will see that the speed of the deflagration front is 
mainly determined by the flame's Arrhenius exponent 



Wf 



u 



knT 



(75) 



that can be large if the energy bias /S.E is small, resulting 
in low Tf. One can express Wf in the form 



1 



fb 



(76) 



The maximal value of Wf compatible with the full- 
burning condition (for = 1) is given by 



W/^max = ^—r, — ~ 6, 



Ah 



fb 



(77) 



where the numerical value corresponds to Mni2 and uses 
Eq. ((74)1 . If Wf exceeds this value, one cannot neglect 

n^_^'^\Tf) in determining the flame temperature. 
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B. The deflagration front 

The deflagration is described by the solution of Eqs. 
P^ . Since in molecular magnets at low temperatures 
both the heat capacity and thermal conductivity strongly 
depend on temperature, it is more convenient to use the 
phonon energy £ as the dynamical variable instead of T, 
see Eq. (fT9| . The speed of the deflagration front can be 
estimated if one rewrites the heat conduction and relax- 
ation equations in a reduced form 



dS ~ dh 

— = V-kV£~ — 

OT OT 

a7 = 



(78) 



which assumes full-burning. The reduced variables are 
defined by 



£ = 



£ 



r 



(79) 



where k/ and Tj are thermal diffusivity and relaxation 
rate at the flame temperature Tf and 



Ir, 



(80) 



describes the width of the deflagration front. The re- 
duced rate is 



T{£) = — = exp 



Wf [1- 



(81) 



where k = k/kj, and T = T /Tf . 

The moving flat deflagration front is a solution of Eqs. 
()78|) that depends on the combined time-like argument 



x/v, 



(82) 



where v is the reduced deflagration speed. In terms of u 
Eqs. ([751) take the form 



d£ 

du 
dfi 

du 



1 d ^d£ 
du du 

= -tn. 



dn 
du 



(83) 



They represent a nonlinear eigenvalue problem with re- 
spect to V. The real deflagration speed v is given by 



(84) 



The advantage of using £ as a dynamical variable is 
that the first of Eqs. (|83|) can be integrated, leading to 



(£+n-l 
du V 



(85) 



Far before and far behind the front £ — const and one 
recovers energy conservation, Eq. (j69p . One can combine 



Eq. (|85p with the second of Eqs. (|83|) to eliminate u. This 
results in 



K{£)f{£)h^ = -i^{£+h-l 
dn 



(86) 



that should be solved with the boundary conditions £ = 
at n = 1 and £ — I at h — 0. 



1. Analytical theory 

Eq. ([85|) allows one to find the spatial variation of tem- 
perature before the front, u < and |m| ^ a, where 
n = 1. For the power-law dependences of Eqs. and 
(fT7| one has k{£) = £-/3/("+i)^ and the solution of Eq. 
([55]) yields the power-law asymptote before the front, 
M < 0, 



£ ^ 



a + 1 



i-u + c) 



-(«+l)//3 



f ^ £l/(« + l) 



.(87) 

where c is the integration constant related to the position 
of the deflagration front. In the realistic case of /3 — 13/3 
the exponent in the expression for T is rather small. That 
is, the heat propagates far ahead of the deflagration front 
due to the divergence of thermal diffisivity at low temper- 
ature. In the case of constant k the temperature before 
the front decreases exponentially as one moves away from 
the front. 



Here the integration constant additive to u was set to 
zero. Now one can find the variation of h before the front 
from the second of Eqs. (|83|) and Eq. (|87|) . to confirm 
that 1 — n is very small. 

Behind the front the T is close to the flame tempera- 
ture, T = £ ~ 1, so that F = 1 and from the second of 
Eqs. ([83)1 one obtains 



Now one can flnd the deviation 6£ 
front from Eq. ([SS]) . Setting k{£) 
a„£ (^£j^p-u _ that yields 



£ = l- 



1 



= 1 - 



1 



(89) 

£—1 behind the 
> 1, one obtains 



(90) 



This relation between £ and h also can be obtained from 
Eq. ([861) with k{£) ^ 1 and f ^ 1 behind the front. 

The speed of the deflagration front can be calculated 
analytically in the high-barrier limit Wf ^1. In this 
case, burning occurs in the region where the temperature 
is already very close to the flame temperature, £ = 1. 
Linearizing the argument of the exponential in Eq. ((ST|) 
on (5£ = £ — 1, one obtains for the relaxation rate 



F^e» 



y = h'fS£ 



(91) 
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with 



vf = Wf 



n-^iAE 



(92) 



The parameter I'f is similar to v of Eq. (|26p. only it is 
defined with respect to the temperature Tf. Using Eqs. 
and ([7T|) . one can simplify J^/ to 



a + 1' 



(93) 



According to Eq. (|77p , the maximal value of j^/ compati- 
ble with the full-burning approximation, is i^/.max ~ 1-5. 
Nevertheless, for simplicity we will consider the case 
v f ^ 1 within the full- burning approximation. In this 
case burning occurs only when the phonon energy is very 
close to its final value, i.e., S£ ^ ^/^f ^ 1- Hence in Eq. 
([85)l one can make a replacement £ ^ 1 and k{£) ^ I'm. 
the burning region. Eq. (j86p then takes the form 



(94) 



dy ~2 -y 

dfi 



It is convenient to consider this equation as an equation 
for n{y). The solution satisfying the boundary conditions 
fi = \ before the front {y — ~oo) and n = behind the 
front {y = 0) reads 



1 



=2/ 



(95) 



It exists if the reduced front speed is given by Vfv'^ — 1 
or 

1 



V 



< 1. 



(96) 



Note that this result is insensitive to the temperature de- 
pendences of the heat capacity and thermal conductivity. 
In real units, one obtains from Eq. ([84]) 



(97) 



where i^/ is given by Eqs. ([5^ or Now 6£ can be 

found from the full system of equations 



dy 
du 



dfi 
du 



(98) 



With account of Vfir^ — 1 and Eq. 195]) the first of these 
equations becomes duy = 1 — e^. The solution is 



y= — In (l + e 
Then from Eq. ([95|l one obtains 

~ -J— 1 ( 

n = — — 

1 + e" 2 V 



u 

1 — tanh — 
2 



(99) 
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FIG. 7: Analytical results for the deflagration front in molec- 
ular magnets in the high-barrier case, v = 0.3 (i.e., Vf = 
l/v^ ~ 11). 




FIG. 8: Numerical results for the reduced deflagration speed 
V vs the flame's Arrhenius exponent Wf for different temper- 
ature dependences Cph oc and k oc . Dashed lines are 
the VK/ > 1 asymptotes v ^ a/(q + 1) /Wf. 



Its accuracy is assured by smallness of w^. In the case of 
K = const, merging Eqs. (l99|) and ([88]) yields 



(102) 



which is the limit of /? ^ in Eq. pUTjl . Note that the 
width of the deflagration front Id defined by Eq. ([80)1 is 
the width of the region where the magnetization changes, 
see Eq. (|100l) . The width of the region where the tem- 
perature changes is Id/ v = ^ HfV f /Y f , according to Eqs. 
([82ll and (fTOS]) . For thermal diffusivity diverging at T = 
for /3 > 0, the region of the temperature variation be- 
comes very broad and its width cannot be defined. 



The solution for £ in the whole range of u can be ob- 
tained by merging Eqs. (|99p and ([57)) . The result is 



£ 



1 



l + -^w2 ln(l + e-") 



(a-H)//3- 



(101) 



2. Numerical results for the deflagration front 



If lyf of Eq. ([92)1 is not large, the problem of the de- 
flagration front cannot be solved analytically. Numeri- 
cal solution uses Eq. ([85]) and the second of Eqs. (j83p . 
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One starts in the region behind the front, with proper 
boundary conditions and arbitrary v, and solves equa- 
tions numerically into the region sufficiently far ahead of 
the front. If the value of v is correct, and only in this 
case, the solution ahead of the front will be n = 1 and 
£ having the asymptotic form of Eq. (l87|) . One finds v 
numerically from one of these conditions, with consistent 
results. The numerically calculated dependences of v on 
Wf are shown in Fig. [S] Surprisingly, the high-barrier 
analytical theory works very well in the whole range of 
Wf. 

To the contrary, the validity of the analytical theory in 
the realistic case of a = 3 and /3 = 13/3 requires rather 
large Wf. One of the reasons for this is that the large- 
Wf approximation in fact requires large i^/ in Eq. ()93|) . 
Nonzero exponents a and (3 have the following effect on 
the deflagration speed. For a > 1 the heat-conduction 
equation is written in terms of £. The decrease of £ ahead 
of the front leads to the decrease of f{£) = that 
enters the relaxation rate F. For a > this temperature 
decrease is less essential than in the case of a = 0. Thus 
the temperature before the front is higher and the front 
moves faster because of faster relaxation. The role of 
/3 > is similar. The heat diffision in the region before 
the front is faster, the temperature before the front is 
higher, and the deflagration speed increases. 



VII. FULL NUMERICAL SOLUTION OF THE 
DEFLAGRATION PROBLEM 

Please, get the full text of the paper here: 
www.lehman.edu/faculty/dgaranin/deflagration.pdf 



VIII. DISCUSSION 

So far experimental work on magnetic deflagration has 
been limited to measurements of the flame speed. The 
proposed theory provides the framework for more de- 
tailed experimental studies suggested in this section. 



A. Deflagration threshold 

According to the theory the magnetic deflagration in a 
crystal of magnetic molecules can be ignited by either in- 
creasing magnetic field or temperature. The simplest sit- 
uation is when the magnetic field and temperature of the 
sample boundary Tq are independent of coordinates. In 
this case the crystal loses stability against formation and 
propagation of the flame (magnetic avalanche) when the 
rate of the spin flip for an individual molecule, T{H, Tq), 
exceeds 



8kiTo)kBT^ 



(103) 



Here fc(To) is coefficient of thermal conductivity at T = 
To, U{H) and AE{H) are field-dependent energy bar- 
rier and energy difference respectively between spin-up 
and spin-down states, I is some characteristic length, 
and n^_i is the initial fraction of molecules available for 
burning. It can be expressed via (negative) initial mag- 
netization Mi and the saturation magnetization Mq as 
n„_i = {Mq — Mi)/ {2Mq). As to the parameter I, it is 
uniquely determined by geometry and is of order of the 
smallest dimension of the crystal. Eq. p03p provides the 
dependence of the critical magnetic field on the temper- 
ature of the sample, or, inversely, the dependence of the 
critical temperature of the sample on the magnetic field, 
see Fig. [21 The deeps in Tq{H) at regularly spaced fields 
are due to the maxima of r(iJ, Tq) at tunneling reso- 
nances. Note that Eq. (|103p contains explicit dependence 
of the deflagration threshold on the initial magnetization 
that should be easy to test in experiment. 

Similar relations have been obtained by us in the pres- 
ence of field and temperature bias, see Sec. IIVI Among 
other things we have demonstrated that the bias sup- 
presses deflagration. The most important outcome of this 
studies is elucidation of the nature of magnetic avalanche. 
Contrary to the initial beliefs, the avalanche does not de- 
velop from a small nucleus of the magnetization reversal 
inside the crystal. It begins as an instability of a smooth 
temperature profile when the spin-flip rate of individual 
molecules (the burning rate) exceeds the rate at which 
the heat flows out of the burning region. The effect is 
exponentially sensitive to the magnetic field and temper- 
ature of the sample. Even a slight fluctuation of H or 
To may take the system deep inside the instability region 
(F 3> Fc), thus, explaining the abrupt and sometimes 
unpredictable nature of the avalanche. It will be inter- 
esting to see if experiments confirms our predictions for 
the deflagration threshold at various initial conditions. 



B. Ignition time 

In this paper we have addressed situations when the 
crystal is instantaneously brought inside the instability 
region. In Sec.|V]we have demonstrated that the ignition 
of the deflagration occurs after a finite time elapses from 
the moment when the instability threshold is crossed. 
When crossing the threshold on field or temperature the 
ignition rate changes from zero below the threshold to 
a some finite value above the threshold, see Figs. ?? - 
??. The deeper one penetrates into the instability re- 
gion the smaller is the ignition time. Most of the experi- 
ments on magnetic avalanches were done in a field-sweep 
mode, when the magnetic field changes at a constant rate, 
H = rt, from a large negative value to a large positive 
value. In this case, the field H'^ at which the avalanche 
occurs should be approximately determined by the equa- 
tion H'^ ^ He + rTig{He), where He is the critical field 
at a temperature Tq plotted in Fig. [2] Since the ignition 
time rapidly falls as H grows above He, it is clear that 
for sufficiently small sweep rates r, the field H'e must be 
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very close to He- This condition is always fulfilled in 
a field-sweep experiment unless a pulse field technique 
is used with a very large r. Thus, a typical field-sweep 
experiment is capable of testing the Hc(Tq) dependence 
plotted in Fig. [2l but not probing the ignition time. To 
measure the ignition time one should apply different tech- 
niques. The trick is to cross the stability threshold by a 
finite step on field or temperature during the time in- 
terval that is small compared to the ignition time. This 
can be achieved by using a small coil with a short time 
constant in addition to the large coil needed to bring the 
system close to the threshold. Alternatively, one can use 
fast heaters to cross the deflagration threshold on the 
temperature of the sample or on the temperature of one 
end of the sample along the lines of Section [Vh. 

Results obtained in one dimension (Sections |V] and 
IVIip show that under symmetric conditions the avalanche 
ignites in the middle of the sample or simultaneously at 
two symmetric regions away from the center. Symme- 
try arguments suggest that this should also be the case 
for any symmetric sample. Meantime experiments done 
under uniform field and temperature conditions often re- 
port the ignition of the deflagration persistently at one 
end of the sample. Explanation of this observation should 
be sought in the inevitable asymmetry of the heat flow 
inside and out of the real crystal. Such an asymmetry oc- 
curs due to the asymmetry of the crystal shape, internal 
inhomogeneity, or as a result of the asymmetric thermal 
insulation from the environment. It should lead to the 
asymmetric temperature profile of a quasi-equilibrium 
state below the deflagration threshold and, thus, asym- 
metric ignition of the deflagration. Our studies elucidate 
the crucial role of the boundary conditions. This should 
be addressed in future experiments by studying, e.g., de- 
flagration in thermally insulated crystals alongside with 
crystals that freely exchange heat with the environment, 
and under asymmetric boundary conditions. 



accuracy to the exponent, the prefactor was estimated by 
order of magnitude. In Section IVII we provided a more 
detailed study of the developed deflagration. Our result 
for the speed of the deflagration front reads 



v{H) 



UkBTfK{Tf)T{H,Tf) 



U{H) 



(104) 



where k, is thermal diffusivity 
Od 



Tf = 



3 ksOo 



1/4 



(105) 



is the flame temperature (the temperature behind the 
front), and Qd is the Debye temperature. Numerical 
exercise with numbers for Mni2 and flelds used in exper- 
iment immediately shows that the above formulas give 
correct estimate of v and Tf and their correct fleld depen- 
dence. Future experiments should show whether these 
formulas provide quantitative description of the devel- 
oped magnetic deflagration. 

An interesting observation that follows from our theory 
is that in the developed magnetic deflagration the width 
of the region where magnetization reverses is different 
from the width of the region inside which the temperature 
decays from Tf to Tq. In fact, the latter region is very 
broad and even difficult to define, see Fig. [7l This is a 
result of the divergence of thermal diffusivity at T ^ 
0, which makes low-temperature magnetic defiagration 
different from chemical defiagration. The latter has a 
well-defined width. Id ^ \J K{Tf) /r{Tf). In the magnetic 
case, however, this formula applies only to the width of 
the region where the magnetization reverses, but not to 
the region where the temperature changes. This may 
explain reported difficulties in local measurements of the 
temperature during the defiagration process. 



C. Velocity and width of the deflagration front 

Measured field and temperature dependences of the 
velocity of the magnetic avalanche are in a reasonably 
good agreement with the concept of deflagration. ^'^^ This 
agreement, however, has only been established with an 
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